C SOURCES:
C    Erich Goll,  Hans-Joachim Werner, and Hermann Stoll.
C    A short-range gradient-corrected densityfunctional in long-range coupled-cluster calculations for rare gas dimers.
C    Physical Chemistry Chemical Physics, 7(23):3917, 2005.

C    Erich Goll, Hans-Joachim Werner, Hermann Stoll, Thierry Leininger, Paola Gori-Giorgi, and Andreas Savin.
C    A short-range gradient-corrected spin density functional in combination with long-range coupled-cluster methods:  Application to alkali-metal rare-gas dimers.
C    Chemical Physics,329(1-3):276-282, oct 2006.



C*****************************************************************************
      pure subroutine ESRX_PBE_GWS_ERF_case_1(rho_a, gamma_aa, mu, Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea
      real*8 :: x0, x1, x2
      Ea = 0.0d0
      x0 = rho_a**1.33333333333333d0
      x1 = 0.217336917462899d0
      x2 = mu**2*rho_a**(-0.666666666666667d0)*x1
      Ea = -1.8610514726982d0*x0*(1.804d0 - 0.804d0/(0.239205142311446d0
     &*gamma_aa* rho_a**(-2.66666666666667d0)*x1*(5.45136177849642d0*x2 
     &- 7.0d0)* exp(-1.4385538026587772d0*x2)/(0.233951463049512d0*mu**4
     &/x0 + 49.1113297378199d0*mu*rho_a**(-0.333333333333333d0) - 49.062
     &2560064678d0*x2 - 81.0d0) + 1.0d0))
      end subroutine


C*****************************************************************************
      pure subroutine ESRX_PBE_GWS_ERF_case_2_1(rho_a, gamma_aa, mu, Ea)
     &
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11
      Ea = 0.0d0
      x0 = rho_a**1.33333333333333d0
      x1 = 4.60115111447049d0
      x2 = 1d0/x1
      x3 = mu**2
      x4 = rho_a**0.666666666666667d0
      x5 = x2*x3/x4
      x6 = rho_a**0.333333333333333d0
      x7 = 1d0/x6
      x8 = 2.14502939711103d0
      x9 = x7/x8
      x10 = mu*x9
      x11 = 0.00844343197019482d0*1d0/rho_a
      Ea = 4.9628039271952d0*x0*(1.804d0 - 0.804d0/(0.239205142311446d0*
     &gamma_aa* rho_a**(-2.66666666666667d0)*x2*(5.45136177849642d0*x5 -
     & 7.0d0)* exp(-1.4385538026587772d0*x5)/(0.233951463049512d0*mu**4/
     &x0 + 49.1113297378199d0*mu*x7 - 49.0622560064678d0*x5 - 81.0d0) + 
     &1.0d0))*(0.275160604074552d0*x10*(mu**3*x11 - mu*(x11*x3 - 0.55032
     &1208149104d0*x9)*exp(-3.3019272488946267d0*x1*x4/x3) - 0.825481812
     &223657d0*x10 + 1.77245385090552d0*erf( 1.8171205928321397d0*x6*x8/
     &mu)) - 0.375d0)
      end subroutine


C*****************************************************************************
      pure subroutine ESRX_PBE_GWS_ERF_case_2_2(rho_a, gamma_aa, mu, Ea)
     &
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11
      Ea = 0.0d0
      x0 = rho_a**1.33333333333333d0
      x1 = 4.60115111447049d0
      x2 = 1d0/x1
      x3 = mu**2
      x4 = 1d0/x3
      x5 = rho_a**0.666666666666667d0
      x6 = x1*x5
      x7 = 2.14502939711103d0
      x8 = rho_a**0.333333333333333d0
      x9 = 1/(x7*x8)
      x10 = mu*x9
      x11 = 0.00844343197019482d0*1d0/rho_a
      Ea = 4.9628039271952d0*x0*(1.804d0 + 0.804d0/(0.239205142311446d0*
     &gamma_aa* rho_a**(-2.66666666666667d0)*x2*x4*(0.213719543995407d0*
     &x0*x4 - 0.183440402716368d0*x6 + 14.4026298885989d0*rho_a**2.0d0/m
     &u**4)* exp(-1.4385538026587772d0*x2*x3/x5) - 1.0d0))*( 0.275160604
     &074552d0*x10*(mu**3*x11 - mu*(x11*x3 - 0.550321208149104d0*x9)*exp
     &(-3.3019272488946267d0*x4*x6) - 0.825481812223657d0*x10 + 1.772453
     &85090552d0*erf( 1.8171205928321397d0*x7*x8/mu)) - 0.375d0)
      end subroutine


C*****************************************************************************
      pure subroutine ESRX_PBE_GWS_ERF_case_2_3(rho_a, gamma_aa, mu, Ea)
     &
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15
      Ea = 0.0d0
      x0 = rho_a**1.33333333333333d0
      x1 = 2.14502939711103d0
      x2 = rho_a**0.333333333333333d0
      x3 = 1/(x1*x2)
      x4 = mu*x3
      x5 = 0.00844343197019482d0*1d0/rho_a
      x6 = erf(1.8171205928321397d0*x1*x2/mu)
      x7 = 4.60115111447049d0
      x8 = mu**2
      x9 = 1d0/x8
      x10 = rho_a**0.666666666666667d0
      x11 = 3.30192724889463d0*x10*x7*x9
      x12 = x8/(x10*x7)
      x13 = exp(x11)
      x14 = 0.0472353356922751d0*mu**4*(x13 - 1.0d0)/x0
      x15 = 3.0d0*x13
      Ea = 4.9628039271952d0*x0*(1.804d0 + 0.804d0/(3.15935190989552d0*g
     &amma_aa* rho_a**(-2.0d0)*x9*(0.00280419844572861d0*x12*(7.0d0*x13 
     &+ 11.0d0 ) - 0.015286700226364d0*x14 + 0.0185185185185185d0)*exp( 
     &-1.4385538026587772d0*x12)/(-0.60570686427738d0*x12*(x15 - 2.0d0 )
     & + 3.90167577847542d0*x13*x4*x6 + 0.183440402716368d0*x14 - x15 ) 
     &- 1.0d0))*(0.275160604074552d0*x4*(mu**3*x5 - mu*( -0.550321208149
     &104d0*x3 + x5*x8)*exp(-x11) - 0.825481812223657d0* x4 + 1.77245385
     &090552d0*x6) - 0.375d0)
      end subroutine


C*****************************************************************************
      pure subroutine ESRX_PBE_GWS_ERF_case_3(rho_a, gamma_aa, mu, Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea
      real*8 :: x0, x1, x2, x3, x4, x5
      Ea = 0.0d0
      x0 = mu**2
      x1 = 1d0/x0
      x2 = rho_a**2.0d0
      x3 = 4.60115111447049d0
      x4 = 1d0/x3
      x5 = rho_a**0.666666666666667d0
      Ea = -3.14159265358979d0*x1*x2*(1.804d0 + 0.804d0/(0.2392051423114
     &46d0* gamma_aa*rho_a**(-2.66666666666667d0)*x1*x4*(0.2137195439954
     &07d0* rho_a**1.33333333333333d0*x1 - 0.183440402716368d0*x3*x5 + 1
     &4.4026298885989d0*x2/mu**4)*exp(-1.4385538026587772d0*x0*x4/x5) - 
     &1.0d0))
      end subroutine


C*****************************************************************************
      pure subroutine D1ESRX_PBE_GWS_ERF_case_1(rho_a, gamma_aa, mu, Ea,
     & d1Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      x0 = 0.682784063255296d0
      x1 = rho_a**1.33333333333333d0
      x2 = 0.826307487110758d0
      x3 = 0.0472353356922751d0
      x4 = 1d0/x1
      x5 = 0.217336917462899d0
      x6 = mu**2
      x7 = rho_a**(-0.666666666666667d0)*x5*x6
      x8 = 1d0/(4.95289087334194d0*mu**4*x3*x4 + 59.4346904801033d0*mu*r
     &ho_a**( -0.333333333333333d0)*x2 - 49.0622560064678d0*x7 - 81.0d0)
     &
      x9 = x8*exp(-1.4385538026587772d0*x7)
      x10 = gamma_aa*x9
      x11 = 5.45136177849642d0*x7 - 7.0d0
      x12 = x11*x5
      x13 = 0.239205142311446d0*rho_a**(-2.66666666666667d0)*x12
      x14 = x10*x13 + 1.0d0
      x15 = 1.804d0 - 0.804d0/x14
      x16 = x14**(-2)
      x17 = rho_a**(-4.33333333333333d0)*x3*x6
      Ea = -2.72568088924821d0*x0*x1*x15
      d1Ea(1) = -x0*(3.63424118566428d0*rho_a**0.333333333333333d0*x15 +
     & 2.19144743495556d0*x1*x10*x16*(mu*x13*x8*(6.60385449778925d0*mu**
     & 3*rho_a**(-2.33333333333333d0)*x3 - 32.7081706709785d0*mu*rho_a**
     & (-1.66666666666667d0)*x5 + 19.8115634933678d0*x2*x4) - 0.63788037
     &9497189d0*rho_a**(-3.66666666666667d0)*x12 + 0.229406311391776d0*x
     &11*x17 - 0.869329180010942d0*x17))
      d1Ea(2) = -0.0777890465506374d0*x11*x16*x4*x9
      end subroutine


C*****************************************************************************
      pure subroutine D1ESRX_PBE_GWS_ERF_case_2_1(rho_a, gamma_aa, mu, E
     &a, d1Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      x0 = 2.14502939711103d0
      x1 = 1d0/x0
      x2 = rho_a**0.333333333333333d0
      x3 = 1d0/x2
      x4 = x1*x3
      x5 = 0.550321208149104d0*x4
      x6 = mu**2
      x7 = 0.00844343197019481d0
      x8 = 1d0/rho_a*x7
      x9 = x6*x8
      x10 = -x5 + x9
      x11 = 4.60115111447049d0
      x12 = rho_a**0.666666666666667d0
      x13 = exp(-3.3019272488946267d0*x11*x12/x6)
      x14 = mu*x13
      x15 = 1d0/mu
      x16 = x0*x15
      x17 = mu**3
      x18 = mu*x4
      x19 = x17*x8 - 0.825481812223657d0*x18 + 1.77245385090552d0*erf( 1
     &.8171205928321397d0*x16*x2)
      x20 = -x10*x14 + x19
      x21 = 0.275160604074552d0*x18
      x22 = x20*x21 - 0.375d0
      x23 = 0.682784063255296d0
      x24 = 0.826307487110758d0
      x25 = 0.0472353356922751d0
      x26 = rho_a**1.33333333333333d0
      x27 = 1d0/x26
      x28 = 1d0/x11
      x29 = x28*x6/x12
      x30 = 1d0/(4.95289087334194d0*mu**4*x25*x27 + 59.4346904801033d0*m
     &u*x24*x3 - 49.0622560064678d0*x29 - 81.0d0)
      x31 = x30*exp(-1.4385538026587772d0*x29)
      x32 = gamma_aa*x31
      x33 = 5.45136177849642d0*x29 - 7.0d0
      x34 = x28*x33
      x35 = 0.239205142311446d0*rho_a**(-2.66666666666667d0)*x34
      x36 = x32*x35 + 1.0d0
      x37 = 1.804d0 - 0.804d0/x36
      x38 = x23*x37
      x39 = 7.26848237132856d0*x26
      x40 = x21*(x14*(x5 - x9) + x19) - 0.375d0
      x41 = x1*x27
      x42 = rho_a**(-2.0d0)*x7
      x43 = x36**(-2)
      x44 = rho_a**(-4.33333333333333d0)*x25*x6
      Ea = x22*x38*x39
      d1Ea(1) = -0.318309886183791d0*mu*x37*x39*(0.0917202013581841d0*x2
     &0*x27 - 0.275160604074552d0*x3*(0.275160604074552d0*mu*x41 + 1.211
     &41372855476d0*rho_a**(-0.666666666666667d0)*x13*x16 + 2.2012848325
     &9642d0*rho_a**(-0.333333333333333d0)*x10*x11*x13*x15 + x14*(-0.183
     &440402716368d0*x41 + x42*x6) - x17*x42)) + 9.69130982843808d0*rho_
     &a**0.333333333333333d0*x38*x40 + 5.84385982654816d0*x23*x26*x32*x4
     &0*x43*(mu*x30*x35*( -32.7081706709785d0*mu*rho_a**(-1.666666666666
     &67d0)*x28 + 6.60385449778925d0*rho_a**(-2.33333333333333d0)*x17*x2
     &5 + 19.8115634933678d0*x24*x27) - 0.637880379497189d0*rho_a**( -3.
     &66666666666667d0)*x34 + 0.229406311391776d0*x33*x44 - 0.8693291800
     &10942d0*x44)
      d1Ea(2) = 0.207437457468366d0*x22*x27*x31*x33*x43
      end subroutine


C*****************************************************************************
      pure subroutine D1ESRX_PBE_GWS_ERF_case_2_2(rho_a, gamma_aa, mu, E
     &a, d1Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      x0 = 0.682784063255296d0
      x1 = 4.60115111447049d0
      x2 = rho_a**0.666666666666667d0
      x3 = x1*x2
      x4 = 0.183440402716368d0*x3
      x5 = rho_a**1.33333333333333d0
      x6 = 21.1705915781930d0
      x7 = mu**2
      x8 = 1d0/x7
      x9 = x6*x8
      x10 = 0.010095114404623d0*x5*x9
      x11 = rho_a**2.0d0
      x12 = 97.4090910340024d0/mu**4
      x13 = 0.147857142857143d0*x11*x12
      x14 = x10 + x13 - x4
      x15 = 1d0/x1
      x16 = x15*x8
      x17 = 0.239205142311446d0*rho_a**(-2.66666666666667d0)*x16
      x18 = exp(-1.4385538026587772d0*x15*x7/x2)
      x19 = gamma_aa*x18
      x20 = x14*x17*x19
      x21 = x20 - 1.0d0
      x22 = 2.14502939711103d0
      x23 = 1d0/x22
      x24 = rho_a**0.333333333333333d0
      x25 = 1d0/x24
      x26 = x23*x25
      x27 = 0.550321208149104d0*x26
      x28 = rho_a**1.0d0
      x29 = 0.00844343197019481d0
      x30 = x29/x28
      x31 = x30*x7
      x32 = -x27 + x31
      x33 = exp(-3.3019272488946267d0*x3*x8)
      x34 = mu*x33
      x35 = 1d0/mu
      x36 = x22*x35
      x37 = mu**3
      x38 = mu*x26
      x39 = x30*x37 - 0.825481812223657d0*x38 + 1.77245385090552d0*erf( 
     &1.8171205928321397d0*x24*x36)
      x40 = -x32*x34 + x39
      x41 = 0.275160604074552d0*x38
      x42 = x40*x41 - 0.375d0
      x43 = 7.26848237132856d0*x5
      x44 = rho_a**0.333333333333333d0
      x45 = -x20 + 1.0d0
      x46 = 1.804d0 - 0.804d0/x45
      x47 = x0*(x41*(x34*(x27 - x31) + x39) - 0.375d0)
      x48 = 1d0/x5
      x49 = x23*x48
      x50 = x29/x11
      x51 = rho_a**(-0.333333333333333d0)*x1
      x52 = -x10 - x13 + x4
      Ea = x0*x42*x43*(1.804d0 + 0.804d0/x21)
      d1Ea(1) = -0.318309886183791d0*mu*x43*x46*(-0.275160604074552d0*x2
     &5*( 0.275160604074552d0*mu*x49 + 1.21141372855476d0*rho_a**( -0.66
     &6666666666667d0)*x33*x36 + 2.20128483259642d0*x32*x33*x35* x51 + x
     &34*(-0.183440402716368d0*x49 + x50*x7) - x37*x50) + 0.091720201358
     &1841d0*x40*x48) + 5.84385982654816d0*x19*x47*x5*( 0.22940631139177
     &6d0*rho_a**(-4.33333333333333d0)*x52/x6 - 0.637880379497189d0*rho_
     &a**(-3.66666666666667d0)*x16*x52 + x17*( -0.295714285714286d0*x12*
     &x28 - 0.0134601525394973d0*x44*x9 + 0.122293601810912d0*x51))/x45*
     &*2 + 9.69130982843808d0*x44*x46*x47
      d1Ea(2) = -0.207437457468366d0*x14*x18*x42*x48*x8/x21**2
      end subroutine


C*****************************************************************************
      pure subroutine D1ESRX_PBE_GWS_ERF_case_2_3(rho_a, gamma_aa, mu, E
     &a, d1Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      x0 = rho_a**1.33333333333333d0
      x1 = 1d0/x0
      x2 = rho_a**0.666666666666667d0
      x3 = 4.60115111447049d0
      x4 = mu**2
      x5 = 1d0/x4
      x6 = x3*x5
      x7 = 3.30192724889463d0*x2*x6
      x8 = exp(x7)
      x9 = 0.0472353356922751d0*mu**4*(x8 - 1.0d0)
      x10 = x1*x9
      x11 = 0.015286700226364d0*x10
      x12 = 7.0d0*x8 + 11.0d0
      x13 = 1d0/x3
      x14 = x13*x4
      x15 = x14/x2
      x16 = 0.00280419844572861d0*x12*x15
      x17 = -x11 + x16 + 0.0185185185185185d0
      x18 = 3.0d0*x8
      x19 = x18 - 2.0d0
      x20 = 2.14502939711103d0
      x21 = 1d0/x20
      x22 = rho_a**0.333333333333333d0
      x23 = 1d0/x22
      x24 = x21*x23
      x25 = mu*x24
      x26 = 1d0/mu
      x27 = x20*x26
      x28 = erf(1.8171205928321397d0*x22*x27)
      x29 = x28*x8
      x30 = 1d0/(0.183440402716368d0*x10 - 0.60570686427738d0*x15*x19 - 
     &x18 + 3.90167577847542d0*x25*x29)
      x31 = rho_a**(-2.0d0)
      x32 = 3.15935190989552d0*x31*x5
      x33 = x30*x32
      x34 = exp(-1.4385538026587772d0*x15)
      x35 = gamma_aa*x34
      x36 = x17*x33*x35
      x37 = x36 - 1.0d0
      x38 = 7.26848237132856d0*x0
      x39 = 0.682784063255296d0
      x40 = 0.550321208149104d0*x24
      x41 = 1d0/rho_a
      x42 = 0.00844343197019481d0
      x43 = x41*x42
      x44 = x4*x43
      x45 = -x40 + x44
      x46 = exp(-x7)
      x47 = mu*x46
      x48 = mu**3
      x49 = -0.825481812223657d0*x25 + 1.77245385090552d0*x28 + x43*x48
      x50 = -x45*x47 + x49
      x51 = 0.275160604074552d0*x25
      x52 = x39*(x50*x51 - 0.375d0)
      x53 = -x36 + 1.0d0
      x54 = 1.804d0 - 0.804d0/x53
      x55 = x39*(x51*(x47*(x40 - x44) + x49) - 0.375d0)
      x56 = x1*x21
      x57 = mu*x56
      x58 = x31*x42
      x59 = rho_a**(-0.666666666666667d0)
      x60 = x27*x59
      x61 = rho_a**(-0.333333333333333d0)
      x62 = x11 - x16 - 0.0185185185185185d0
      x63 = x41*x8
      x64 = rho_a**(-1.66666666666667d0)*x14
      x65 = x64*x8
      x66 = rho_a**(-2.33333333333333d0)*x9
      Ea = x38*x52*(1.804d0 + 0.804d0/x37)
      d1Ea(1) = -0.318309886183791d0*mu*x38*x54*(0.0917202013581841d0*x1
     &*x50 - 0.275160604074552d0*x23*(2.20128483259642d0*x26*x3*x45*x46*
     &x61 + 1.21141372855476d0*x46*x60 + x47*(x4*x58 - 0.183440402716368
     &d0* x56) - x48*x58 + 0.275160604074552d0*x57)) + 9.69130982843808d
     &0* rho_a**0.333333333333333d0*x54*x55 + 5.84385982654816d0*x0*x30*
     & x35*x55*(3.02993180261165d0*rho_a**(-3.66666666666667d0)*x13*x62 
     &- 6.31870381979104d0*rho_a**(-3.0d0)*x5*x62 + x32*( 0.001869465630
     &48574d0*x12*x64 - 0.0432098765432099d0*x63 + 0.0336503813487433d0*
     &x65 - 0.0203822669684853d0*x66) + x33*x62*( -0.40380457618492d0*x1
     &9*x64 + 1.30055859282514d0*x29*x57 - 8.58869971286677d0*x29*x60 - 
     &2.66666666666667d0*x41 + 6.60385449778925d0*x6*x61*x8 + 4.0d0*x63 
     &- 0.40380457618492d0*x65 + 0.244587203621824d0*x66))/x53**2
      d1Ea(2) = -18.4628097041666d0*x17*x30*x34*x5*x52*x59/x37**2
      end subroutine


C*****************************************************************************
      pure subroutine D1ESRX_PBE_GWS_ERF_case_3(rho_a, gamma_aa, mu, Ea,
     & d1Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      x0 = rho_a**2.0d0
      x1 = 4.60115111447049d0
      x2 = rho_a**0.666666666666667d0
      x3 = 0.183440402716368d0*x1*x2
      x4 = mu**2
      x5 = 1d0/x4
      x6 = 21.1705915781930d0
      x7 = x5*x6
      x8 = 0.010095114404623d0*rho_a**1.33333333333333d0*x7
      x9 = mu**(-4)
      x10 = 97.4090910340024d0*x9
      x11 = 0.147857142857143d0*x0*x10
      x12 = x11 - x3 + x8
      x13 = 1d0/x1
      x14 = x13*x5
      x15 = 0.239205142311446d0*rho_a**(-2.66666666666667d0)*x14
      x16 = exp(-1.4385538026587772d0*x13*x4/x2)
      x17 = gamma_aa*x16
      x18 = x12*x15*x17
      x19 = x18 - 1.0d0
      x20 = 3.14159265358979d0*x5
      x21 = rho_a**1.0d0
      x22 = -x18 + 1.0d0
      x23 = -x11 + x3 - x8
      Ea = -1.0d0*x0*x20*(1.804d0 + 0.804d0/x19)
      d1Ea(1) = -x20*(0.804d0*x0*x17*(0.229406311391776d0*rho_a**(-4.333
     &33333333333d0)* x23/x6 - 0.637880379497189d0*rho_a**(-3.6666666666
     &6667d0)*x14*x23 + x15*(0.122293601810912d0*rho_a**(-0.333333333333
     &333d0)*x1 - 0.0134601525394973d0*rho_a**0.333333333333333d0*x7 - 0
     &.295714285714286d0*x10*x21))/x22**2 + 2.0d0*x21*(1.804d0 - 0.804d0
     &/x22))
      d1Ea(2) = 0.131313669051252d0*rho_a**(-0.666666666666667d0)*x12*x1
     &6*x9/x19**2
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRX_PBE_GWS_ERF_case_1(rho_a, gamma_aa, mu, Ea,
     & d1Ea, d2Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4), d2Ea(10)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      d2Ea(:) = 0.0d0
      x0 = 0.217336917462899d0
      x1 = mu**2
      x2 = x0*x1
      x3 = rho_a**(-0.666666666666667d0)*x2
      x4 = exp(-1.4385538026587772d0*x3)
      x5 = gamma_aa*x4
      x6 = 0.826307487110758d0
      x7 = 0.0472353356922751d0
      x8 = mu**4
      x9 = rho_a**1.33333333333333d0
      x10 = 1d0/x9
      x11 = 59.4346904801033d0*mu*rho_a**(-0.333333333333333d0)*x6 + 4.9
     &5289087334194d0*x10*x7*x8 - 49.0622560064678d0*x3 - 81.0d0
      x12 = 1d0/x11
      x13 = rho_a**(-2.66666666666667d0)
      x14 = 5.45136177849642d0*x3 - 7.0d0
      x15 = x0*x14
      x16 = x13*x15
      x17 = 0.239205142311446d0*x16
      x18 = x12*x17
      x19 = x18*x5 + 1.0d0
      x20 = 1d0/x19
      x21 = -0.804d0*x20 + 1.804d0
      x22 = 0.682784063255296d0
      x23 = x22*x9
      x24 = rho_a**0.333333333333333d0
      x25 = rho_a**(-4.33333333333333d0)
      x26 = x1*x7
      x27 = x25*x26
      x28 = rho_a**(-3.66666666666667d0)*x15
      x29 = x14*x27
      x30 = x10*x6
      x31 = mu*x0
      x32 = rho_a**(-1.66666666666667d0)*x31
      x33 = rho_a**(-2.33333333333333d0)
      x34 = mu**3*x7
      x35 = x33*x34
      x36 = 19.8115634933678d0*x30 - 32.7081706709785d0*x32 + 6.60385449
     &778925d0*x35
      x37 = x12*x36
      x38 = mu*x37
      x39 = x17*x38 - 0.869329180010942d0*x27 - 0.637880379497189d0*x28 
     &+ 0.229406311391776d0*x29
      x40 = x39*x4
      x41 = 2.19144743495556d0*x40
      x42 = x12/x19**2
      x43 = gamma_aa*x42
      x44 = 0.148394183600699d0*x14
      x45 = x4*x44
      x46 = 1.27576075899438d0*x28
      x47 = 2.19144743495556d0*x9
      x48 = x11**(-2)
      x49 = exp(-2.8771076053175544d0*x3)
      x50 = x48*x49/x19**3
      x51 = 0.0102659822546843d0*rho_a**(-6.0d0)*x8
      x52 = rho_a**(-5.33333333333333d0)*x26
      x53 = x25*x34*x37
      Ea = -2.72568088924821d0*x21*x23
      d1Ea(1) = -x22*(3.63424118566428d0*x21*x24 + x41*x43*x9)
      d1Ea(2) = -0.524205495546598d0*x10*x42*x45
      d2Ea(1) = x22*(gamma_aa**2*x39*x47*x50*(0.478410284622892d0*x16*x3
     &8 - 1.73865836002188d0*x27 + 0.458812622783553d0*x29 - x46) - 1.21
     &141372855476d0*rho_a**(-0.666666666666667d0)*x21 - 5.8438598265481
     &6d0*x24*x40*x43 - x42*x47*x5*(mu*x18*( -15.4089938281749d0*rho_a**
     &(-3.33333333333333d0)*x34 + 54.5136177849642d0*x13*x31 - 26.415417
     &991157d0*x33*x6) + 2.33889472482303d0*rho_a**(-4.66666666666667d0)
     &*x15 + 0.239205142311446d0*x13*x14*x2*x36*x48*(39.6231269867355d0*
     &x30 - 65.416341341957d0*x32 + 13.2077089955785d0*x35) + 0.22000888
     &1071042d0*x14*x51 - 1.60584417974243d0*x14*x52 + 0.458812622783553
     &d0*x14*x53 - x38*x46 - 1.66743573022264d0*x51 + 6.08530426007659d0
     &*x52 - 1.73865836002188d0*x53))
      d2Ea(2) = x42*(1.0484109910932d0*gamma_aa*x10*x12*x20*x39*x44*x49 
     &- x23*x41 - 0.698940660728797d0*x33*x45)
      d2Ea(3) = 0.00808821074878272d0*rho_a**(-4.0d0)*x14**2*x50
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRX_PBE_GWS_ERF_case_2_1(rho_a, gamma_aa, mu, E
     &a, d1Ea, d2Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4), d2Ea(10)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      d2Ea(:) = 0.0d0
      x0 = 2.14502939711103d0
      x1 = 1d0/x0
      x2 = rho_a**0.333333333333333d0
      x3 = 1d0/x2
      x4 = x1*x3
      x5 = 0.550321208149104d0*x4
      x6 = mu**2
      x7 = 1d0/rho_a
      x8 = 9.86960440108936d0
      x9 = 1d0/x8
      x10 = 0.0833333333333333d0*x9
      x11 = x10*x7
      x12 = x11*x6
      x13 = x12 - x5
      x14 = 4.60115111447049d0
      x15 = rho_a**0.666666666666667d0
      x16 = exp(-3.3019272488946267d0*x14*x15/x6)
      x17 = mu*x16
      x18 = 1d0/mu
      x19 = x0*x18
      x20 = mu**3
      x21 = mu*x4
      x22 = x11*x20 - 0.825481812223657d0*x21 + 1.77245385090552d0*erf( 
     &1.8171205928321397d0*x19*x2)
      x23 = -x13*x17 + x22
      x24 = 0.275160604074552d0*x21
      x25 = x23*x24 - 0.375d0
      x26 = 0.682784063255296d0
      x27 = 1d0/x14
      x28 = x27*x6
      x29 = x28/x15
      x30 = exp(-1.4385538026587772d0*x29)
      x31 = gamma_aa*x30
      x32 = 0.826307487110758d0
      x33 = 21.1705915781930d0
      x34 = 1d0/x33
      x35 = mu**4
      x36 = rho_a**1.33333333333333d0
      x37 = 1d0/x36
      x38 = 59.4346904801033d0*mu*x3*x32 - 49.0622560064678d0*x29 + 4.95
     &289087334194d0*x34*x35*x37 - 81.0d0
      x39 = 1d0/x38
      x40 = rho_a**(-2.66666666666667d0)
      x41 = 5.45136177849642d0*x29 - 7.0d0
      x42 = x27*x41
      x43 = x40*x42
      x44 = 0.239205142311446d0*x43
      x45 = x39*x44
      x46 = x31*x45 + 1.0d0
      x47 = 1d0/x46
      x48 = -0.804d0*x47 + 1.804d0
      x49 = x26*x48
      x50 = x25*x49
      x51 = 7.26848237132856d0*x36
      x52 = rho_a**0.333333333333333d0
      x53 = x24*(x17*(-x12 + x5) + x22) - 0.375d0
      x54 = 0.318309886183791d0
      x55 = mu*x37
      x56 = rho_a**(-2.0d0)*x10
      x57 = rho_a**(-0.666666666666667d0)
      x58 = x16*x19
      x59 = 0.183440402716368d0*x37
      x60 = -x1*x59 + x56*x6
      x61 = rho_a**(-0.333333333333333d0)*x14*x16*x18
      x62 = 0.275160604074552d0*x1*x55 + 2.20128483259642d0*x13*x61 + x1
     &7*x60 - x20* x56 + 1.21141372855476d0*x57*x58
      x63 = 0.275160604074552d0*x3
      x64 = 0.0917202013581841d0*x23*x37 - x62*x63
      x65 = mu*x54*x64
      x66 = x48*x51
      x67 = x39/x46**2
      x68 = x30*x53
      x69 = rho_a**(-4.33333333333333d0)
      x70 = x34*x6
      x71 = x69*x70
      x72 = rho_a**(-3.66666666666667d0)*x42
      x73 = x41*x71
      x74 = x32*x37
      x75 = mu*x27
      x76 = rho_a**(-1.66666666666667d0)*x75
      x77 = rho_a**(-2.33333333333333d0)
      x78 = x20*x34
      x79 = x77*x78
      x80 = 19.8115634933678d0*x74 - 32.7081706709785d0*x76 + 6.60385449
     &778925d0*x79
      x81 = x39*x80
      x82 = mu*x81
      x83 = x44*x82 - 0.869329180010942d0*x71 - 0.637880379497189d0*x72 
     &+ 0.229406311391776d0*x73
      x84 = x26*x83
      x85 = 5.84385982654816d0*x36
      x86 = x84*x85
      x87 = x68*x86
      x88 = 1.39788132145759d0*x30
      x89 = x25*x67
      x90 = 0.148394183600699d0*x41
      x91 = x37*x90
      x92 = x31*x89
      x93 = 1.27576075899438d0*x72
      x94 = x38**(-2)
      x95 = exp(-2.8771076053175544d0*x29)
      x96 = x25*x94*x95/x46**3
      x97 = x1*x77
      x98 = 0.166666666666667d0*rho_a**(-3.0d0)*x9
      x99 = x16/x20
      x100 = 0.0102659822546843d0*rho_a**(-6.0d0)*x35
      x101 = rho_a**(-5.33333333333333d0)*x70
      x102 = x69*x78*x81
      Ea = x50*x51
      d1Ea(1) = gamma_aa*x67*x87 + 9.69130982843808d0*x49*x52*x53 - x65*
     &x66
      d1Ea(2) = x88*x89*x91
      d2Ea(1) = -gamma_aa**2*x86*x96*(0.478410284622892d0*x43*x82 - 1.73
     &865836002188d0* x71 + 0.458812622783553d0*x73 - x93) + mu*x54*x66*
     &( 0.122293601810912d0*x23*x77 - x59*x62 + x63*(-0.366880805432736d
     &0 *mu*x97 - 0.80760915236984d0*rho_a**(-1.66666666666667d0)*x58 - 
     &0.733761610865473d0*rho_a**(-1.33333333333333d0)*x13*x14*x16*x18 -
     & 4.84565491421904d0*x13*x33*x57*x99 - x17*(x6*x98 - 0.244587203621
     &824d0*x97) + x20*x98 - 4.40256966519284d0*x60*x61 - 2.666666666666
     &67d0*x7*x8*x99)) + x26*x85*x92*(-mu*x45*( 15.4089938281749d0*rho_a
     &**(-3.33333333333333d0)*x78 + 26.415417991157d0*x32*x77 - 54.51361
     &77849642d0*x40*x75) + 2.33889472482303d0*rho_a**(-4.66666666666667
     &d0)*x42 + 0.220008881071042d0*x100*x41 - 1.66743573022264d0*x100 -
     & 1.60584417974243d0*x101*x41 + 6.08530426007659d0*x101 + 0.4588126
     &22783553d0*x102*x41 - 1.73865836002188d0*x102 + 0.239205142311446d
     &0*x28*x40*x41*x80*x94*(39.6231269867355d0*x74 - 65.416341341957d0*
     &x76 + 13.2077089955785d0*x79) - x82*x93) - 11.6877196530963d0*x31*
     &x36*x65*x67*x83 - 19.3826196568762d0*x48* x52*x65 + 3.230436609479
     &36d0*x50*x57 + 15.5836262041284d0*x52*x84 *x92
      d2Ea(2) = x67*(-2.79576264291519d0*gamma_aa*x39*x47*x53*x83*x91*x9
     &5 - 0.0691804894611514d0*x41*x55*x64*x88 + 1.86384176194346d0*x68*
     &x77 *x90 + x87)
      d2Ea(3) = -0.0215685619967539d0*rho_a**(-4.0d0)*x41**2*x96
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRX_PBE_GWS_ERF_case_2_2(rho_a, gamma_aa, mu, E
     &a, d1Ea, d2Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4), d2Ea(10)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104
     &, x105, x106
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      d2Ea(:) = 0.0d0
      x0 = rho_a**1.33333333333333d0
      x1 = 7.26848237132856d0*x0
      x2 = mu**2
      x3 = 1d0/x2
      x4 = 4.60115111447049d0
      x5 = 1d0/x4
      x6 = x3*x5
      x7 = rho_a**(-2.66666666666667d0)*x6
      x8 = 0.239205142311446d0*x7
      x9 = rho_a**0.666666666666667d0
      x10 = x4*x9
      x11 = 0.183440402716368d0*x10
      x12 = 21.1705915781930d0
      x13 = x12*x3
      x14 = 0.010095114404623d0*x0*x13
      x15 = rho_a**2.0d0
      x16 = 97.4090910340024d0
      x17 = mu**(-4)
      x18 = x16*x17
      x19 = 0.147857142857143d0*x15*x18
      x20 = -x11 + x14 + x19
      x21 = x2*x5/x9
      x22 = exp(-1.4385538026587772d0*x21)
      x23 = x20*x22
      x24 = gamma_aa*x23*x8
      x25 = x24 - 1.0d0
      x26 = 1.804d0 + 0.804d0/x25
      x27 = 0.682784063255296d0
      x28 = 2.14502939711103d0
      x29 = 1d0/x28
      x30 = rho_a**0.333333333333333d0
      x31 = 1d0/x30
      x32 = x29*x31
      x33 = 0.550321208149104d0*x32
      x34 = rho_a**1.0d0
      x35 = 1d0/x34
      x36 = 9.86960440108936d0
      x37 = 1d0/x36
      x38 = 0.0833333333333333d0*x37
      x39 = x35*x38
      x40 = x2*x39
      x41 = -x33 + x40
      x42 = exp(-3.3019272488946267d0*x10*x3)
      x43 = mu*x42
      x44 = 1d0/mu
      x45 = x28*x44
      x46 = mu**3
      x47 = mu*x32
      x48 = x39*x46 - 0.825481812223657d0*x47 + 1.77245385090552d0*erf( 
     &1.8171205928321397d0*x30*x45)
      x49 = -x41*x43 + x48
      x50 = 0.275160604074552d0*x47
      x51 = x49*x50 - 0.375d0
      x52 = x27*x51
      x53 = x26*x52
      x54 = rho_a**0.333333333333333d0
      x55 = -x24 + 1.0d0
      x56 = 1d0/x55
      x57 = -0.804d0*x56 + 1.804d0
      x58 = x50*(x43*(x33 - x40) + x48) - 0.375d0
      x59 = x27*x58
      x60 = 1d0/x0
      x61 = x38/x15
      x62 = rho_a**(-0.666666666666667d0)
      x63 = x42*x45
      x64 = 0.183440402716368d0*x60
      x65 = x2*x61 - x29*x64
      x66 = rho_a**(-0.333333333333333d0)*x4
      x67 = x42*x44*x66
      x68 = 0.275160604074552d0*mu*x29*x60 + 2.20128483259642d0*x41*x67 
     &+ x43*x65 - x46*x61 + 1.21141372855476d0*x62*x63
      x69 = 0.275160604074552d0*x31
      x70 = 0.0917202013581841d0*x49*x60 - x68*x69
      x71 = 0.318309886183791d0
      x72 = mu*x1*x71
      x73 = x55**(-2)
      x74 = x11 - x14 - x19
      x75 = 1d0/x12
      x76 = rho_a**(-4.33333333333333d0)*x75
      x77 = 0.229406311391776d0*x76
      x78 = rho_a**(-3.66666666666667d0)*x6
      x79 = 0.637880379497189d0*x78
      x80 = 0.122293601810912d0*x66
      x81 = 0.0134601525394973d0*x13*x54
      x82 = 0.295714285714286d0*x18
      x83 = x34*x82
      x84 = x74*x77 - x74*x79 + x8*(x80 - x81 - x83)
      x85 = 5.84385982654816d0*x0
      x86 = x22*x85
      x87 = x59*x84*x86
      x88 = x25**(-2)
      x89 = 1.39788132145759d0*x23*x60
      x90 = 0.148394183600699d0*x3
      x91 = mu*x70*x71
      x92 = -x80 + x81 + x83
      x93 = -x20*x77 + x20*x79 - x8*x92
      x94 = x52*x93
      x95 = gamma_aa*x88
      x96 = x22*x95
      x97 = 0.458812622783553d0*x76
      x98 = rho_a**(-1.33333333333333d0)*x4
      x99 = 1.27576075899438d0*x78
      x100 = exp(-2.8771076053175544d0*x21)
      x101 = x100/x25**3
      x102 = rho_a**(-2.33333333333333d0)
      x103 = x102*x29
      x104 = 0.166666666666667d0*rho_a**(-3.0d0)*x37
      x105 = x42/x46
      x106 = x58*x90
      Ea = x1*x53
      d1Ea(1) = gamma_aa*x73*x87 + 9.69130982843808d0*x54*x57*x59 - x57*
     &x70*x72
      d1Ea(2) = -x51*x88*x89*x90
      d2Ea(1) = gamma_aa**2*x101*x85*x94*(-x20*x97 + x20*x99 - 0.4784102
     &84622892d0*x7* x92) - 11.6877196530963d0*x0*x91*x93*x96 - 19.38261
     &96568762d0*x26 *x54*x91 + x26*x72*(0.122293601810912d0*x102*x49 - 
     &x64*x68 + x69* (-0.366880805432736d0*mu*x103 - 0.80760915236984d0*
     &rho_a**( -1.66666666666667d0)*x63 + x104*x46 - 4.84565491421904d0*
     &x105*x12 *x41*x62 - 2.66666666666667d0*x105*x35*x36 - 0.7337616108
     &65473d0* x41*x42*x44*x98 - x43*(-0.244587203621824d0*x103 + x104*x
     &2) - 4.40256966519284d0*x65*x67)) + x52*x86*x95*(-0.22000888107104
     &2d0* rho_a**(-6.0d0)*x2*x20/x16 + 1.60584417974243d0*rho_a**( -5.3
     &3333333333333d0)*x20*x75 - 2.33889472482303d0*rho_a**( -4.66666666
     &666667d0)*x20*x6 - x8*(0.00448671751316577d0*x13*x62 + x82 + 0.040
     &7645339369707d0*x98) - x92*x97 + x92*x99) + 3.23043660947936d0*x53
     &*x62 + 15.5836262041284d0*x54*x94*x96
      d2Ea(2) = x73*(2.79576264291519d0*gamma_aa*x100*x106*x20*x56*x60*x
     &84 - 1.86384176194346d0*x102*x106*x23 + 0.0691804894611514d0*x44*x
     &70* x89 + x87)
      d2Ea(3) = 0.0215685619967539d0*rho_a**(-4.0d0)*x101*x17*x20**2*x51
     &
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRX_PBE_GWS_ERF_case_2_3(rho_a, gamma_aa, mu, E
     &a, d1Ea, d2Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4), d2Ea(10)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104
     &, x105, x106, x107, x108, x109, x110, x111, x112, x113, x114, x115
     &, x116, x117, x118, x119, x120, x121, x122, x123, x124, x125, x126
     &, x127, x128, x129, x130, x131, x132, x133, x134, x135, x136, x137
     &, x138, x139, x140, x141, x142
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      d2Ea(:) = 0.0d0
      x0 = rho_a**1.33333333333333d0
      x1 = 7.26848237132856d0*x0
      x2 = rho_a**0.666666666666667d0
      x3 = 4.60115111447049d0
      x4 = mu**2
      x5 = 1d0/x4
      x6 = x3*x5
      x7 = 3.30192724889463d0*x2*x6
      x8 = exp(x7)
      x9 = 3.0d0*x8
      x10 = 1d0/x0
      x11 = 21.1705915781930d0
      x12 = 1d0/x11
      x13 = mu**4
      x14 = x12*x13*(x8 - 1.0d0)
      x15 = x10*x14
      x16 = x9 - 2.0d0
      x17 = 1d0/x3
      x18 = x17*x4
      x19 = x18/x2
      x20 = 2.14502939711103d0
      x21 = 1d0/x20
      x22 = rho_a**0.333333333333333d0
      x23 = 1d0/x22
      x24 = x21*x23
      x25 = mu*x24
      x26 = 1d0/mu
      x27 = x20*x26
      x28 = erf(1.8171205928321397d0*x22*x27)
      x29 = x28*x8
      x30 = 0.183440402716368d0*x15 - 0.60570686427738d0*x16*x19 + 3.901
     &67577847542d0*x25*x29 - x9
      x31 = 1d0/x30
      x32 = 0.015286700226364d0*x15
      x33 = 7.0d0*x8 + 11.0d0
      x34 = 0.00280419844572861d0*x19*x33
      x35 = -x32 + x34 + 0.0185185185185185d0
      x36 = rho_a**(-2.0d0)
      x37 = x36*x5
      x38 = 3.15935190989552d0*x37
      x39 = x35*x38
      x40 = x31*x39
      x41 = exp(-1.4385538026587772d0*x19)
      x42 = gamma_aa*x41
      x43 = x40*x42
      x44 = x43 - 1.0d0
      x45 = 1.804d0 + 0.804d0/x44
      x46 = 0.682784063255296d0
      x47 = 0.550321208149104d0*x24
      x48 = 1d0/rho_a
      x49 = 9.86960440108936d0
      x50 = 1d0/x49
      x51 = 0.0833333333333333d0*x50
      x52 = x48*x51
      x53 = x4*x52
      x54 = -x47 + x53
      x55 = exp(-x7)
      x56 = mu*x55
      x57 = mu**3
      x58 = -0.825481812223657d0*x25 + 1.77245385090552d0*x28 + x52*x57
      x59 = -x54*x56 + x58
      x60 = 0.275160604074552d0*x25
      x61 = x46*(x59*x60 - 0.375d0)
      x62 = x45*x61
      x63 = rho_a**0.333333333333333d0
      x64 = -x43 + 1.0d0
      x65 = 1d0/x64
      x66 = -0.804d0*x65 + 1.804d0
      x67 = x46*(x60*(x56*(x47 - x53) + x58) - 0.375d0)
      x68 = 0.318309886183791d0
      x69 = mu*x10*x21
      x70 = x36*x51
      x71 = rho_a**(-0.666666666666667d0)
      x72 = x27*x55
      x73 = 0.183440402716368d0*x10
      x74 = -x21*x73 + x4*x70
      x75 = rho_a**(-0.333333333333333d0)
      x76 = x26*x3*x55*x75
      x77 = 2.20128483259642d0*x54*x76 + x56*x74 - x57*x70 + 0.275160604
     &074552d0*x69 + 1.21141372855476d0*x71*x72
      x78 = 0.275160604074552d0*x23
      x79 = x68*(0.0917202013581841d0*x10*x59 - x77*x78)
      x80 = mu*x1
      x81 = x31/x64**2
      x82 = x32 - x34 - 0.0185185185185185d0
      x83 = rho_a**(-3.66666666666667d0)*x17
      x84 = 3.02993180261165d0*x83
      x85 = rho_a**(-3.0d0)
      x86 = x5*x85
      x87 = 6.31870381979104d0*x86
      x88 = x48*x8
      x89 = 0.0432098765432099d0*x88
      x90 = rho_a**(-1.66666666666667d0)
      x91 = x18*x90
      x92 = x8*x91
      x93 = 0.0336503813487433d0*x92
      x94 = rho_a**(-2.33333333333333d0)
      x95 = x14*x94
      x96 = 0.0203822669684853d0*x95
      x97 = 0.00186946563048574d0*x33*x91
      x98 = 0.564189583547756d0
      x99 = x48*x98
      x100 = x6*x8
      x101 = 6.60385449778925d0*x100
      x102 = x16*x91
      x103 = x29*x69
      x104 = x71*x8
      x105 = x104*x27*x28
      x106 = x101*x75 - 0.40380457618492d0*x102 + 1.30055859282514d0*x10
     &3 - 8.58869971286677d0*x105 + 4.0d0*x88 - 0.40380457618492d0*x92 +
     & 0.244587203621824d0*x95 - 4.72654360241471d0*x99
      x107 = x106*x31
      x108 = x107*x38
      x109 = x108*x82 + x38*(-x89 + x93 - x96 + x97) + x82*x84 - x82*x87
     &
      x110 = x41*x67
      x111 = 5.84385982654816d0*x0
      x112 = x109*x110*x111
      x113 = x35*x5
      x114 = 18.4628097041666d0*x41*x71
      x115 = x31/x44**2
      x116 = x115*x61
      x117 = mu*x79
      x118 = x21*x94
      x119 = mu*x118
      x120 = 0.166666666666667d0*x50*x85
      x121 = rho_a**(-1.66666666666667d0)
      x122 = 1d0/x57
      x123 = x122*x55
      x124 = rho_a**(-1.33333333333333d0)
      x125 = x89 - x93 + x96 - x97
      x126 = -x108*x35 - x125*x38 - x35*x84 + x35*x87
      x127 = x116*x42
      x128 = 6.0598636052233d0*x83
      x129 = x128*x35
      x130 = 12.6374076395821d0*x86
      x131 = x130*x35
      x132 = 6.31870381979104d0*x37
      x133 = x125*x132
      x134 = x30**(-2)
      x135 = exp(-2.8771076053175544d0*x19)
      x136 = x134*x135*x61/x44**3
      x137 = x36*x8
      x138 = rho_a**(-2.66666666666667d0)*x18
      x139 = x138*x8
      x140 = rho_a**(-3.33333333333333d0)*x14
      x141 = 1d0/x13
      x142 = x27*x29
      Ea = x1*x62
      d1Ea(1) = gamma_aa*x112*x81 + 9.69130982843808d0*x63*x66*x67 - x66
     &*x79*x80
      d1Ea(2) = -x113*x114*x116
      d2Ea(1) = gamma_aa**2*x111*x126*x136*(-x107*x132*x35 - x129 + x131
     & - x133) - 11.6877196530963d0*x0*x115*x117*x126*x42 + x111*x127*( 
     &-2.90581327762917d0*rho_a**(-5.33333333333333d0)*x12*x35*x4 + 17.1
     &696135481327d0*rho_a**(-4.66666666666667d0)*x17*x35 - 18.956111459
     &3731d0*rho_a**(-4.0d0)*x113 - x106*x134*x39*( 13.2077089955785d0*x
     &100*x75 - 0.80760915236984d0*x102 + 2.60111718565028d0*x103 - 17.1
     &773994257335d0*x105 + 8.0d0*x88 - 0.80760915236984d0*x92 + 0.48917
     &4407243648d0*x95 - 9.45308720482942d0*x99) - x107*x129 + x107*x131
     & - x107*x133 - x125*x128 + x125*x130 - x38*(0.0951172458529316d0*x
     &100*x124 - 0.146090534979424d0*x137 + 0.00311577605080957d0*x138*x
     &33 + 0.10095114404623d0*x139 - 0.0475586229264658d0*x140) - x40*(x
     &101* x124 + 14.5369647426571d0*x104*x11*x141 - 1.73407812376686d0*
     &x119 *x29 + 5.72579980857785d0*x121*x142 - 18.9061744096588d0*x122
     &*x28 *x49*x88 - 27.0091844287146d0*x124*x5 - 7.55555555555556d0*x1
     &37 + 0.673007626974866d0*x138*x16 + 1.21141372855476d0*x139 - 0.57
     &070347511759d0*x140 + 2.86289990428892d0*x142*x90 + 6.302058136552
     &94d0*x36*x98)) - 19.3826196568762d0*x117*x45*x63 + 15.583626204128
     &4d0*x126*x127*x63 + x45*x68*x80*( 0.122293601810912d0*x59*x94 - x7
     &3*x77 + x78*(-4.84565491421904d0* x11*x123*x54*x71 - 0.36688080543
     &2736d0*x119 + x120*x57 - 0.80760915236984d0*x121*x72 - 2.666666666
     &66667d0*x123*x48*x49 - 0.733761610865473d0*x124*x26*x3*x54*x55 - x
     &56*( -0.244587203621824d0*x118 + x120*x4) - 4.40256966519284d0*x74
     &*x76 )) + 3.23043660947936d0*x62*x71
      d2Ea(2) = x81*(36.9256194083333d0*gamma_aa*x109*x113*x135*x31*x65*
     &x67*x71 - 24.6170796055555d0*x110*x113*x121 + x112 + x114*x26*x35*
     &x79)
      d2Ea(3) = 116.661026201793d0*rho_a**(-2.66666666666667d0)*x136*x14
     &1*x35**2
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRX_PBE_GWS_ERF_case_3(rho_a, gamma_aa, mu, Ea,
     & d1Ea, d2Ea)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 21, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_a, gamma_aa, mu
      real*8, intent(out) :: Ea, d1Ea(4), d2Ea(10)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50
      Ea = 0.0d0
      d1Ea(:) = 0.0d0
      d2Ea(:) = 0.0d0
      x0 = rho_a**2.0d0
      x1 = mu**2
      x2 = 1d0/x1
      x3 = 4.60115111447049d0
      x4 = 1d0/x3
      x5 = x2*x4
      x6 = rho_a**(-2.66666666666667d0)
      x7 = 0.239205142311446d0*x6
      x8 = x5*x7
      x9 = rho_a**0.666666666666667d0
      x10 = x1*x4/x9
      x11 = exp(-1.4385538026587772d0*x10)
      x12 = 0.183440402716368d0*x3*x9
      x13 = 21.1705915781930d0
      x14 = x13*x2
      x15 = 0.010095114404623d0*rho_a**1.33333333333333d0*x14
      x16 = 97.4090910340024d0
      x17 = mu**(-4)
      x18 = x16*x17
      x19 = 0.147857142857143d0*x0*x18
      x20 = -x12 + x15 + x19
      x21 = x11*x20
      x22 = gamma_aa*x21*x8
      x23 = x22 - 1.0d0
      x24 = 3.14159265358979d0
      x25 = x2*x24
      x26 = rho_a**1.0d0
      x27 = -x22 + 1.0d0
      x28 = 1d0/x27
      x29 = x12 - x15 - x19
      x30 = 1d0/x13
      x31 = rho_a**(-4.33333333333333d0)*x30
      x32 = x29*x31
      x33 = rho_a**(-3.66666666666667d0)
      x34 = x29*x5
      x35 = x33*x34
      x36 = 0.295714285714286d0*x18
      x37 = 0.122293601810912d0*rho_a**(-0.333333333333333d0)*x3 - 0.013
     &4601525394973d0*rho_a**0.333333333333333d0*x14 - x26*x36
      x38 = x37*x5
      x39 = 0.229406311391776d0*x32 - 0.637880379497189d0*x35 + x38*x7
      x40 = x11*x39
      x41 = 0.804d0*x0
      x42 = x40*x41
      x43 = x27**(-2)
      x44 = gamma_aa*x43
      x45 = rho_a**(-0.666666666666667d0)
      x46 = 0.682784063255296d0
      x47 = x21*x46
      x48 = exp(-2.8771076053175544d0*x10)
      x49 = x39*x48
      x50 = 0.384641868836805d0*x2
      Ea = -1.0d0*x0*x25*(1.804d0 + 0.804d0/x23)
      d1Ea(1) = -x25*(2.0d0*x26*(-0.804d0*x28 + 1.804d0) + x42*x44)
      d1Ea(2) = 0.192320934418403d0*x17*x45*x47/x23**2
      d2Ea(1) = x25*(gamma_aa**2*x41*x49*(0.458812622783553d0*x32 - 1.27
     &576075899438d0* x35 + 0.478410284622892d0*x38*x6)/x27**3 - x11*x41
     &*x44*( 0.220008881071042d0*rho_a**(-6.0d0)*x1*x29/x16 - 1.60584417
     &974243d0*rho_a**(-5.33333333333333d0)*x29*x30 + 2.33889472482303d0
     &*rho_a**(-4.66666666666667d0)*x34 + 0.458812622783553d0*x31*x37 - 
     &1.27576075899438d0*x33*x38 - x8*( 0.0407645339369707d0*rho_a**(-1.
     &33333333333333d0)*x3 + 0.00448671751316577d0*rho_a**(-0.6666666666
     &66667d0)*x14 + x36)) - 3.216d0*x26*x40*x44 + 1.608d0*x28 - 3.608d0
     &)
      d2Ea(2) = x2*x43*(-gamma_aa*x20*x28*x45*x46*x49*x50 + rho_a**(-1.6
     &6666666666667d0) *x47*x50 - x24*x42)
      d2Ea(3) = -0.0136534984882616d0*rho_a**(-3.33333333333333d0)*x20**
     &2*x48/(mu**6*x23 **3)
      end subroutine
